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ABSTRACT 

Context. In dense parts of interstellar clouds (> 10^ cm“^), dust and gas are expected to be in thermal equilibrium, being coupled via 
collisions. However, previous studies have shown that in the presence of intense radiation fields, the temperatures of the dust and gas 
may remain decoupled even at higher densities. 

Aims. The objective of this work is to study in detail the temperatures of dust and gas in the photon-dominated region S 140, especially 
around the deeply embedded infrared sources IRS 1-3 and at the ionization front. 

Methods. We derive the dust temperature and column density by combining Herschel PACS continuum observations with SOFIA 
observations at 37 pm and SCUBA data at 450 pm. We model these observations using simple greybody fits and the DUSTY radiative 
transfer code. For the gas analysis we use RADEX to model the CO 1-0, CO 2-1, *^CO 1-0 and C**0 1-0 emission lines mapped with 
the IRAM-30m telescope over a 4' field. Around IRS 1-3, we use HIFI observations of single-points and cuts in CO 9-8, *^CO 10-9 
and C*®0 9-8 to constrain the amount of warm gas, using the best fitting dust model derived with DUSTY as input to the non-local 
radiative transfer model RATRAN. The velocity information in the lines allows us to separate the quiescent component from outflows 
when deriving the gas temperature and column density. 

Results. We find that the gas temperature around the infrared sources varies between ~35 and ~55 K. In contrast to expectation, the 
gas is systematically warmer than the dust by ~5-15 K despite the high gas density. In addition we observe an increase of the gas 
temperature from 30-35 K in the surrounding up to 40-45 K towards the ionization front, most likely due to the UV radiation from 
the external star. Furthermore, detailed models of the temperature structure close to IRS 1 which take the known density gradient into 
account show that the gas is warmer and/or denser than what we model. Finally, modelling of the dust emission from the sub-mm 
peak SMM 1 constrains its luminosity to a few x 10^ Lq. 

Conclusions. We conclude that the gas heating in the S 140 region is very efficient even at high densities. The most likely explanation 
is deep UV penetration from the embedded sources in a clumpy medium and/or oblique shocks. 

Key words. ISM: individual (S 140), ISM: kinematics and dynamics, ISM: molecules, stars: formation 


1. Introduction 

Stars are born inside dense molecular structures in the inter¬ 
stellar medium (ISM) which consist of gas and dust. The dust 
grains in clouds associated with star-forming regions absorb the 
short-wavelength radiation from the central stars, heat up, and 
then re-emit radiation at far-infrared (FIR) and sub-millimeter 
wavelengths. The thermal radiation from dust is optically thin at 
these wavelengths and thus is a good tracer of physical param¬ 
eters such as temperature, density and gas mass of the clouds. 
Dust is efficiently heated through near IR radiation due to the 
broadband absorption properties of the dust grains. On the other 
hand, the gas temperature Tg^^, is mainly governed by heating 
processes driven by UV radiation (photoelectric heating, H 2 ex¬ 
citation, dissociation). 


* Based on Herschel observations. Herschel is an ESA space obser¬ 
vatory with science instruments provided by European-led Principal 
Investigator consortia and with important participation from NASA. 


In the well-shielded centers of massive cores, the primary 
effects for the thermal balance are cosmic ray heating, molecu¬ 
lar line cooling and collisional heating or cooling due to dust- 
gas collisions (Goldsmith & Langer 1978; Tielens 2005; Draine 
2011 ). 

At densities above ~10^'^, the dust and gas are expected to 
be collisionally coupled and they are characterized by the same 
temperature (Goldsmith 2001). At lower densities, the rate of 
collisions between dust and gas decreases and the cooling of 
the gas via fine-structure and molecular rotational transitions 
becomes dominant (e.g. CO, [C II] and [OI]) (Sternberg & 
Dalgarno 1995; Kaufman et al. 1999; Meijerink & Spaans 2005). 

Excitation of the rotational levels of molecules like CO (low 
dipole moment) observed at (sub)mm wavelengths is provided 
by collisions with H 2 . If these collisions are frequent enough 
to exceed the spontaneous decay rate, the levels will get into 
thermal equilibrium with H 2 . The critical densities n^r of CO 
and isotopologues for the 1=1-0 and 2-1 transitions at 50 K 
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are ~2xl0^ and ~2xl0'^ respectively At this point 
their excitation temperature will approach the gas kinetic tem¬ 
perature and thus the study of these lines is ideal for gas kinetic 
temperature estimates. The rarer isotopologues of CO, such as 
and C'^O, can be used in order to probe regions of high 
column density where the lines of the abundant isotopologues 
become optically thick. The critical densities n^r of CO and 
isotopologues for the J=9-8 and 10-9 transitions at 50 K are 
-1.2x10^ cm“^ and ~1.5xl0® cm“^ respectively. 

S 140 is a well studied HII region that lies at the south¬ 
west edge of the molecular cloud L 1204 at a distance of 746 pc 
(Hirota et al. 2008). In previous studies (e.g. Dedes et al. 2010; 
Ikeda & Kitamura 2011) a larger distance of 910 pc (Crampton 
& Fisher 1974) has been adopted. This edge is illuminated by the 
BOV star HD211880, creating a visible HII region and a Photon- 
Dominated Region (PDR) at the ionization front (IF). The cloud 
also hosts a cluster of embedded high-mass young stellar objects 
(YSOs), IRS 1 to 3 (Evans et al. 1989) at a projected distance of 
~ 75" northeast of the IF. A large number of low-mass stars are 
also forming, making S 140 an ideal laboratory for studying star 
formation of various masses and the differences caused by irra¬ 
diation from external (BOV star) and internal sources (IRS 1 to 
3). 

The analysis of physical models around IRS 1 has been per¬ 
formed by a number of previous authors including Harvey et al. 
(1978), Giirtler et al. (1991), Minchin et al. (1993), van der Tak 
et al. (2000), Mueller et al. (2002), de Wit et al. (2009), Maud 
& Hoare (2013). Poelman & Spaans (2006) have reported the 
dumpiness towards S 140 giving a value of n^j ~ 10“* cm“^ 
for the interclump medium and ~4xl0^ for the clump gas. 
Previous models lacked spatial information at many wavelengths 
and many of them were based on the assumption that the gas 
temperature is coupled to the dust as a result of the high densi¬ 
ties (e.g. Poelman & Spaans 2006). 

In this study we have analyzed the temperature and column 
density around the S 140 high-mass star-forming region of em¬ 
bedded young stellar objects by combining Herschel PACS and 
HIFI data with ground-based mapping observations of mm- 
wave lines using the IRAM 30 m telescope. The new PACS con¬ 
tinuum data and IRAM spectroscopic maps provide us with the 
highest angular resolution spatial information available cover¬ 
ing an area that includes the YSOs in both datasets. The IF was 
covered by our larger IRAM maps (4'x4') but not by the PACS 
observations, since the latter covers a smaller area (45"x45"). 
Furthermore, the data give us the opportunity to study inde¬ 
pendently the gas and dust temperatures around FIRS 1-3 and 
check whether the assumption of well coupled dust-gas is valid 
throughout this cloud. 

In the following sections, we present the dust and gas mod¬ 
eling, and the resulting parameters. We present the observations 
in §2, the observational results in §3, the gas - dust analysis and 
the comparison between gas and dust in §4, the effect of density 
gradients using more advanced dust and gas models in §5 and a 
discussion of our results and conclusions in §6. 

2. Observations and data reduction 

2.1. PACS data 

To investigate the temperature and column-density structure of 
the dust we have analyzed a single footprint of PACS (Poglitsch 

* The values were calculated using Einstein coefficients and colli- 
sional rates (Eq. (2)) from Yang et al. (2010). 


et al. 2010) with its 5x5 spatial array of ~ 9" pixels over a range 
of wavelengths from ~ 70 fxm to ~ 200 /rm (obsid: 1342222256). 
The data reduction was performed using HIPE 9.1. 

Since the PACS spectrometer continuum data cover only the 
peak of the spectral energy distribution (SED) for S 140, we have 
supplemented our analysis with shorter and longer wavelength 
images. We used the SOEIA/EORCAST images at 11.1, 31.5, 
and 37 /rm discussed by Harvey et al. (2012) and the 24.5 /urn 
Subaru/COMICS data of de Wit et al. (2009). Since all of these 
observations were obtained with resolution better than that of the 
PACS/Spec images, we used them both at their native resolution 
and re-convolved to the PACS/Spec 187 fim resolution. These 
raw and convolved images were also re-sampled to a 1" grid like 
the PACS/Spec images. Longward of the PACS wavelengths the 
highest resolution image available is that in the JCMT SCUBA 
archive^ at 450 /rm. We have used this image without further 
processing in our analysis below. We have also examined the 
SCUBA 850 /rm image which is qualitatively similar but with 
lower spatial resolution. 


2.2. IRAM data 

The 30m observations were conducted during 4 hours on 
September 24th, 2011. The Eight Mixer Receiver (EMIR) was 
used during the observations. CO 7 = 1 —> 0 and isotopologues 
(~ 110-115 GHz, 3 mm) and CO 7 = 2 ^ 1 (230.54 GHz, 
1 mm) map data were gathered in two linear polarizations. Our 
4'x4' maps were centered close to S 140-IRS 1 position at RA 
= 22:19:18.30, Dec = 63:18:54.2 (J2000) (Eigure 2). Eor this 
work, we selected the maps of CO 7 = 1 ^0(115.27 GHz), 7 
= 2^ 1 (230.54 GHz) together with the isotopologues '^CO 7 
= 1^0 (110.20 GHz) and C'*0 7=1^ 0 (109.78 GHz). At 
all the observed lines, the ETS and WILMA were connected in 
parallel as backends. In the 3 mm set up we collected data from 
3 bands of ~ 4 GHz each including the isotopic lines of ^^CO 7 
= 1 —» 0 and C'*0 7 = 1 —> 0 in the CO setup. The beam sizes 
are ~ 21" at 3 mm and ~ 11" at 1 mm. 

The sky opacity measured by the taumeter at 225 GHz was 
stable at T 225 = -0.26. Pointing and focus were checked about 
every hour on quasars and on Uranus, and were stable within 
2” and 0.2 mm. The spectral line calibration was checked by 
pointed observations on IRS 1 and IP. The on-the-fly map¬ 
ping observations were done in perpendicular directions to avoid 
scanning artefacts. The typical system temperatures were 96- 
215 K at 3 mm, and 175 - 190 K at 1 mm. The resulting rms 
falls between 0.152-0.705 K at 3 mm and ~ 0.525 K at 1 mm 
for a velocity resolution of ~ 0.5 km s“^ and ~0.25 km s“' re¬ 
spectively. 

The intensities were converted from antenna temperature 
units to a scale of main-beam temperature (T^b), dividing by 
the main-beam efficiencies rjmb i^eff/Peff) of 0.82 and 0.64 at 
3 mm and 1 mm respectively. Velocities are given with respect 
to the LSR. Baselines of orders 0 up to 4 have been subtracted. 
Data reduction and analysis were performed using the GILDAS 
software ^. 


^ m96bu47199708210025 from http://www.cadc-ccda.hia-iha.nrc- 
cnrc. gc. ca/en/j cmt/ 

^ GILDAS is a radio-astronomy software developed by IRAM. See 
http://www.iram.fr/TRAMFR/GILDAS/ 
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2.3. HI FI data 

In order to test the reliability of our method for the gas analy¬ 
sis, we use single pointing observations (spectra) (Figure 6) and 
cuts from Herschel-HIFI (Figure 2, top), in addition to the 4’ x4’ 
maps from IRAM-30m. Single-pointing, frequency-switching 
spectra in bands la, lb, 2b, 3a, 4a, 4b, 6b and 7b and OTF 
maps were analyzed. The single pointing data were taken to¬ 
wards S 140-IRS 1 at RA = 22:19:18.21, Dec = 63:18:46.9 
(J2000) and towards a position in the IF at RA = 22:19:11.53, 
Dec = 63:17:46.9 (J2000). 

For this analysis, we used only the observations of CO J - 
9 —> 8 (1036.9 GHz), together with the isotopologues, '^CO J - 
10^9 (1101.35 GHz), C'^O 7 = 9^ 8 (987.56 GHz), towards 
IRS 1 and IF (obsids: 1342195050, 1342196426, 1342219209, 
1342195049, 1342201741). In addition we used the cuts of CO 
9^8 and '^CO 7 = 10 ^ 9 (obsids:1342201741, 1342201806). 
The observed line parameters in units of antenna temperatures, 
have been corrected in units of main beam temperatures using 
the main beam efficiency of 0.74 (987.56 GHz, 1036.9 GHz, 
1101.35 GHz) (Roelfsema et al. 2012), while the beam size at 
these frequencies is ~ 21". 


3. Observational results 

3.1. Dust continuum observations 

For this analysis we selected our 45"x45"maps in a number 
of continuum wavelengths with no obvious line emission or 
with line emission that was easily masked. The particular wave¬ 
lengths chosen were centered on the wavelengths for which con¬ 
volution kernels have been developed by Aniano et al. (2011), 
which enables us to convolve all maps to a common resolution. 
These wavelengths PACS/Spec data are 73, 75, 84, 94, 110, 125, 
136, 145, 150, 168, and 187 /um. We extracted 3 % bandwidth 
continuum fluxes centered on each of these wavelengths while 
masking out obvious line emission. For comparison with other 
wavelength data as well as for comparison with modeling re¬ 
sults, we have re-sampled all these images to a 1 "spatial grid 
(Figure 1). For the final comparison the model images are con¬ 
volved to a telescope PSF and typically a 9" square pixel to sim¬ 
ulate the PACS/Spec observations. 

Figure 1 shows the image at 73//m. The 37/zm (SOFIA) 
emission is overlaid using black contours showing the positions 
of IRS 1, 2, and 3. As shown in the image, the emission from 
IRS 1 dominates at this wavelength. 

3.2. The spatiai and veiocity distribution of CO emission 

Figure 2 shows the maps of the peak intensities as observed in 
their original angular resolution at their peak velocities of CO 
1-0, CO 2-1, '^CO 1-0 and C'^O 1-0. The lines of the main 
isotope show two clearly disjunct peaks, one at the densest clus¬ 
ter between the YSOs IRS 1 and IRS 3, and one marking the ex¬ 
ternal interface of the cloud illuminated by HD211880. As '^CO 
and C'^O trace the column density structure, the interface peak 
is shifted into the cloud for these lines, merging with the peak 
around the embedded cluster. 

Finally, while CO 1-0 and CO 2-1 peak towards IRS 1-3 
(ARA: 10.70", ADec: -4.0"), the '^CO 1-0 and C^*^0 1-0 peak 
closer to IRS 2 (ARA: -11.20", ADec: 18.30"). Ci*^0 1-0 shows 
a second peak southeast of IRS 1. This reflects a column density 
effect which is also revealed later in the column density map (see 
Figure 11). 



22"19"’22= 22"19"’18’ 22"l9™14' 


Fig. 1 PACS/Spec continuum image at 73 pm: 5x5 spatial array of 
9" pixels re-sampled to 1" grid with contours of 37pm (SOFIA) emis¬ 
sion overlaid showing the positions of IRS 1, 2, and 3. 


Figure 3 shows channel maps of CO 1-0 to illustrate the ve¬ 
locity structure. We find a clear velocity offset between the gas 
around the cluster and the interface and a much narrower veloc¬ 
ity distribution of the interface gas. Figures 4-5 shows the line 
profiles of CO and isotopologues as observed towards IRS 1 and 
IF. The peak velocity of the emission changes with position in 
the map. The IF peaks close to the velocity of the source which is 
-7.1 km s ', while IRS 1 peaks at-3.6 km s *. The differ¬ 
ence may be caused by outflows driven by the infrared sources. 
For the above reasons we chose to model the peak intensities of 
the lines in their peak velocities. 

3.3. Lineprofiies 

The line profiles as shown in Figures 4 and 5, are stronger to¬ 
wards IRS 1 and weaker towards the IF (Table 1). At the indi¬ 
vidual positions, the lines from CO are broader than the lines 
from the other isotopologues of CO (Figure 4 and Figure 6). 

Table 1 presents the line parameters as measured towards the 
IRS 1 and IF positions, applying a single component Gaussian 
fit, including both HIFI (Figure 6) and IRAM (Figures 4 and 5) 
observations. The uncertainties quoted for the peak intensities 
correspond to the observational RMS and the uncertainties for 
the FWHM and i/lsr are from the Gaussian fits to the peak posi¬ 
tion. The width of the lines varies throughout the cloud, showing 
broader profiles towards the center of the map, where the three 
infrared sources are located and narrow towards the ionization 
front (Figure 7). Being interested in the quiescent gas and thus 
the narrow component of the lines, we use observed and com¬ 
puted peak intensities in our calculations, as the narrow com¬ 
ponent dominates the total peak intensity of the lines (>50%). 
With this method we limit the effects from outflow activities on 
our calculations but we cannot totally exclude them. The outflow 
contamination of the line is stronger towards the central sources 
where the integrated contribution of the broad and narrow com¬ 
ponents is comparable and weaker in positions away from the 
sources where the narrow component provides 70-85% of the 
total peak intensity. We restricted the peak intensity analysis to a 
velocity window of the width of the FWHM of the lines around 
the '^CO 1-0 line. In this way we always model the same quies¬ 
cent component. 
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Table 1 Line parameters of CO and isotopologues towards the IRS 1 and IF positions - IRAM & HlFl data. 



1 IRS 1 

1 IF (-45" ,-60") 

Molecule Transition Rest Frequency 

E,,;, 

ulsr 

FWHM 

T 

mb 

UlSR 

FWHM 

T 

mb 



(MHz) 

(K) 

(km s“') 

(km s“') 

(K) 

(km s“') 

(km s“') 

(K) 

CO 

1-0 

115271.20 

5.53 

-5.07 ± 0.08 

10.82 ±0.20 31.91 ±0.55 

-7.32 ±0.19 

2.66 ± 0.07 

5.80 ±0.08 

CO 

2-1 

230538.00 

16.60 

-5.0 ±0.1 

11.40 ±0.22 

31.69 ±0.65 

-7.77 ± 0.01 

3.20 ± 0.02 

7.62 ± 0.06 

CO 

9-8 

1036912.39 

248.88 

-6.60 ± 0.02 

8.32 ± 0.07 

27.62 ± 0.37 

-7.80±0.07 

2.50±0.20 

1.20±0.20 

'^co 

1-0 

110201.35 

5.29 

-6.87 ± 0.03 

4.09 ± 0.08 

16.52 ± 0.34 

-8.12 ±0.37 

1.34 ± 0.70 

0.84 ± 0.26 

'^co 

10-9 

1101349.60 

290.79 

-6.98 ± 0.03 

3.80 ± 0.07 

8.56 ± 0.03 



<3 RMS 

c'^o 

1-0 

109782.17 

5.27 

-5.71 ±0.14 

2.37 ± 0.39 

2.88 ±0.37 



<3 RMS 

c^'o 

9-8 

987560.20 

237.02 

-6.94 ± 0.06 

3.26 ± 0.04 

1.31 ±0.06 



<3 RMS 


4. Physical conditions 

4.1. RADEX fitting 

4.1.1. Method 

We use the non-LTE radiative transfer program RADEX (van 
der Tak et al. 2007) to compare the observed line intensity ra¬ 
tios with a grid of models for deriving kinetic temperatures, 
gas densities and column densities. As model input we use the 
molecular data from the LAMDA database (Schoier et al. 2005) 
and CO collisional rate coefficients from Yang et al. (2010). 
RADEX predicts line intensities of a molecule for a given set 
of parameters: kinetic temperature, column density, H 2 density, 
background temperature and line width. 

In the density range relevant for S 140 the CO 1-0/2-1 ratio 
is a good tracer for kinetic temperatures. This is demonstrated in 
Eigure 8 showing the CO 1-0/2-1 line ratios for a CO column 
density of 10'^ cm“^. Eor typical gas temperatures, we see that 
these ratios are insensitive to H 2 density for values > 10“* ^ cm“^, 
well above the critical density of both transitions. At this point 
we should also point out that at low kinetic temperatures this 
ratio is also not a good indicator of kinetic temperature at the ac¬ 
curacy needed to distinguish between 30 K and 60 K. In this tem¬ 
perature range CO 1-0/2-1 ~ 1.0 + 0.1, one would have to know 
reliably the main beam temperature to within 10 % to trust the 
solution to the kinetic temperature to within + 15 K. Given the 
uncertainties in antenna and receiver calibrations, beam sizes, 
sidelobes, etc. this level of accuracy is difficult to achieve. One 
main source of uncertainty is the coupling efficiency of the beam 
to the source structure. Eor sources larger than the main beam, 
the 30 m errorbeams (and sidelobes) contribute to the coupling. 
At 230 GHz, the main beam efficiency of the 30 m is high, 60 %, 
and the three 30 m errorbeams contribute 26% to the total power 
received by the beam pattern. In addition, assuming a reasonable 
source size of about 1 ', only the main beam and 1st errorbeam 
contribute. As the 1st errorbeam contains only 4% of the total 
beam power, the correction factor for the antenna temperatures, 
Eeff/Beff decreases only little, from 92/59=1.56 to 92/63=1.46. 
However, it is not the purpose of our paper to de-convolve the 
channel maps for the errorbeams or to provide a full accounting 
of the error budget contributing to T„,h. 

4.1.2. Assumptions 

Eor the CO 1-0, CO 2-1, ^^CO and C^^O peak intensities 
(~ 3x RMS) we performed a;^f^ minimization to fit the gas ki¬ 
netic temperatures and column densities assuming that all lines 
arise from the same gas. Eor more accurate modeling the CO 
2-1 map was convolved to a lower angular resolution in or¬ 
der to be consistent with the other lines (~21"). The func¬ 


tion was computed as the quadratic sum of the differences 
between the observed and the synthetic line intensities for a 
range of kinetic temperatures (10 K<Tk,„<200 K) and column 
densities (10^^cm“^<Nco<10*^cm~^), values that are consistent 
with the expected ones for such region (e.g. van der Tak et al. 
2000; Spaans & van Dishoeck 1997; Poelman & Spaans 2006; 
Hiittemeister et al. 1993). 

An additional free parameter in RADEX is the tempera¬ 
ture of the background radiation field that may pump the line 
transitions. We adopted the value of 2.73 K for all our calcu¬ 
lations since the cosmic background radiation (CMB) peaks at 
1.871 mm and it is the prominent component at millimeter wave¬ 
lengths. We used a fixed line width value of 3.5 km s ' that 
approximates the value that we have measured throughout the 
cloud for the narrow component. Einally, in order to compute 
the intensities of the isotopic lines we assume fixed isotopic ra¬ 
tios of '^CG/'^CG = 60 and C“’G/C'*'G = 560 (Wannier 1980; 
Wilson & Rood 1994). 

4.1.3. Gas density 

Previous studies have shown that the S 140 region contains gas 
with n/fj ~ 10^ cm“^. Poelman & Spaans (2006) reports dumpi¬ 
ness in the region giving a value of n^j ~ 10"' cm“^ for the inter¬ 
clump medium and ~4xl0^ for the clump gas. 

Pig. 8 shows that it is impossible to derive the gas density 
from the low-7 lines of CO (2-1, 1-0) for densities of n^j ~ 
10^ cm“^ or above. A reliable determination of the gas density 
was thus only possible for those parts of the map where the HIPI 
cuts provided additional high-7 line data. Here, we performed a 
three-parameter RADEX analysis, fitting the CO 1-0, CO 2-1, 
CO 9-8, '^CO 1-0, '^CO 10-9 and C'*0 1-0 lines convolved 
to the same resolution (21"), to determine the 11 ^^. The result 
is shown in Pig 9. This fit proves that, at least along the cuts, 
we find everywhere gas densities well above n/^j = 10^ cm“^, 
where the CO 1-O/CO 2-1 line ratio can be directly used as a 
temperature measure. In this way we can be sure to avoid the 
regimes where the ratio is sensitive to the H 2 density. 

To assess the reliability of the results, we perform two sep¬ 
arate RADEX fits along the cuts where the 9-8 and 10-9 tran¬ 
sitions of CO and the isotopologues were observed by HIPI. In 
the first fitting approach we fix n^j to 10^ cm“^ and fit only the 
low-7 CO lines observed by IRAM (blue curve in Eigure 10) 
and all the lines (red curve in Eigure 10). The second fit uses 
all the lines and treats the gas density as a free parameter in the 
range between 7xl0^cm“^ and 7xl0^cm“^ as in Pig. 9. The best 
fit of the latter occurred for n//^ >10^ cm“^. Adding the higher-7 
lines drives the fit to systematically higher kinetic temperatures 
by about ~5-15 K (Eigure 10). The opacities of the lines are 


4 



Koumpia et al: Temperatures of Dust and Gas in S 140. 




100 50 0 -50 -100 


ARA ["1 
1-0 



100 50 0 -50 -100 

ARA ["1 


Fig. 2 Spatial distribution of CO 1-0, CO 2-1, ‘^CO 1-0 & C‘*0 1-0 
peak intensities (FIIFI cut & PACS overplotted). The IRS 1-3 positions 
were taken from Evans et al. (1989) & SMM 1 from Maud et al. (2013). 
The center is atRA = 22:19:18.30, Dec = 63:18:54.2 (J2000). 


found to be in the following ranges: CO 1-0: 16^5, CO 2-1; 
52-154,00 9-8: 3.6-65, '^CO 1-0:0.3-0.8, '^CO 10-9:0.017- 
0.38, C'^O 1-0: 0.025-0.09. As an attempt to test how much the 
optically thin '^CO 10-9 line influences the solution we re-ran 
our calculations applying a double weight for this line. The re- 


Fig. 3 Channel maps of CO 1-0 obtained with IRAM, with the central 
velocities given in km s“*. The levels of contours are 45 K (darker) till 
10 K (lighter) using a step of 5 K. 



^LSR (l^^/sec) 

Fig. 4 Line profiles of CO 1-0 (hlack line), CO 2-1 (grey line), '^CO 
1-0 (blue line) and C'*0 1-0 (red line) towards IRS 1 as observed with 
IRAM (offset position: 0", -4.0"). Wings are visible at blueshifted ve¬ 
locities. The dashed vertical line represent the velocity of the source 
(-7.1 km s“*). 


suiting kinetic temperatures changed by K that is inside the 
range of the reported errors. 

We ran the same kind of analysis towards the two positions 
with most lines observed"^ (see Table 1). Towards IRS 1 the full 
analysis provides a gas kinetic temperature of ~ 60^® K and a 
column density of ~ 3^2Xl0'®cm“^ while the same procedure 


For the low-i CO data we extracted the points from the IRAM maps 
closest to HIFI pointing observations with no more than 3" offset. 
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Fig. 5 Line profiles of CO 1-0 (black line), CO 2-1 (grey line), '^CO 
1-0 (blue line) and C**0 1-0 (red line) towards the IF as observed with 
IRAM (offset position: -42.7", -68.0"). The dashed vertical line repre¬ 
sent the velocity of the source (-7.1 km s“'). 


Fig. 7 Line width as observed throughout the cloud. The FWHMs were 
derived from a 2nd moment map of the line emission. Broader pro¬ 
files appear towards the center (>10 km s“*), where the three in¬ 
frared sources are located and narrower towards the ionization front 
(<5 km s“'). The FWHM that characterizes the gas that surrounds the 
protostars varies between ~5-10 km s“'. 



V (km/sec) 

LSR ^ ' 


Fig. 6 Line profiles of CO 9-8 (black line), '^CO 10-9 (blue line) 
and C'*0 9-8 (red line) towards IRS 1 as observed with HlFl (RA = 
22:19:18.21, Dec = 63:18:46.9). The dashed vertical line represent the 
velocity of the source (-7.1 km s“'). 


CO 1-0 / CO 2-1 



n(H,) [cm"’] 

Fig. 8 Line ratios of CO 1-O/CO 2-1 as a function of kinetic temperature 
and H 2 density for column density N^o = lO'* cm“^. 


using only the IRAM data results in a kinetic temperature of 
~ 49^3 K and a column density of ~ 1 3^ jXlO'^cm^^. Towards 
the IF position, the full dataset indicates a kinetic temperature of 
~ 40^2 K and column density of ~ 1.0;[j jXl0'*cm“^, while the 
IRAM data only result in a kinetic temperature of ~ 35^g K and 
column density of ~ 4.4^2^xlO'^cm”^. The gas temperatures 
obtained when using only the low-7 lines from IRAM underes¬ 
timate Tgflj, providing a lower limit of gas temperatures (Yildiz 
et al. 2013). The major and minor axis of the resultedcontours 
were used as the error bars of the two parameters. 

For a more accurate determination of the H 2 density, tracers 
such as CS 2-1/3-2 and HCO^ 1-0/3-2 are more reliable (van 
der Tak et al. 2007), but were not observed. Snell et al. (1984) 
derived a density of 5xl0^cm“^ towards IRS 1-3 using four CS 
transitions which is in agreement with the value we derive. In 
addition Goldsmith & Danger (1999) performed a LTE popula¬ 
tion diagram analysis towards IRS 1 using the same CS dataset 
and derived a kinetic temperature between 30 K-50 K which is 


lower than our value probably due to the effect of their larger 
beam size. 


4.1.4. Kinetic temperature & column density distribution 

The resulting maps of kinetic temperatures and CO column den¬ 
sities are shown in Figure 11. We find a kinetic temperature of 
~ 55 K toward the center and ~ 45 K toward the ionization front, 
while the rest of the cloud is characterized by lower temperatures 
(25^0 K) (Figure 11). The column density toward the center 
was found to be the highest, with a value of ~ 6.5xl0'*cm“^, 
while toward the IF was found to be ~ l.lxlO'^cm^^. The low¬ 
est column density that was determined throughout the cloud is 
~ 1.4xl0'®cm“^. The gas temperature map obviously reflects 
the different heating contributions, by the embedded cluster of 
the three high-mass YSOs (IRS 1 to 3) in the center and from 
the external BOV star HD211880 towards the IF. 
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Fig. 9 Profiles of gas temperatures (bottom), column densities (middle) 
and n^j densities (top) along the HlFl cut, as resulting from the 3-free 
parameter fitting process to the 30 m and FIlFl data. The resulting n//^ 
densities are always > 10^ cm“^ throughout the HlFl cut. 


4.2. Dust Analysis 

Prior to any radiative transfer modeling of the dust in this re¬ 
gion (Sect.5) we here estimate the rough temperatures and op¬ 
tical depths relevant to the dust distribution. For this we used 
the PACS/Spec continuum data from 73 - 187 /tm which cov¬ 
ers well the peak of the SED over the entire region mapped with 
PACS (Figure 12). To analyze our dust continuum observations, 
we used two subsets of the PACS and SOFIA images. To com¬ 
pute a luminosity map, we used the 11 - 187 fim images all as 
re-convolved to the 187 fim resolution (13"). At each F'pixel of 
the computed image we integrated the flux density from 11.1- 
187 jum and included a linear extrapolation to zero flux level be¬ 
yond 187 fj.m, which typically resulting in only a few % addition 
to the total. We believe the absolute uncertainties in this map 
are of order 15% based on the absolute calibration uncertain¬ 
ties of all the input data. The relative uncertainties are probably 
much less, + 5%, since there is very little luminosity shortward 
of 10 jum and longward of 200 jum. 

We computed color-temperature/optical depth values at each 
spatial point on the 1" grid by fitting a blackbody function mod¬ 
ified by a T dust emissivity variation. The peak dust optical 
depths at the shortest wavelength are likely several xO.l, but 
still less than unity, so we did not include any effects from opti¬ 
cal thickness in these estimates. A comparison of the luminosity, 
color temperature, and optical depth maps and several images is 
shown in Figure 12. These images show a number of important 
qualitative properties of the dust emission. First, there is a clear 
and relatively smooth change in source morphology from that at 
73 jum, which is similar to the 10-37 fim images, to the mor¬ 
phology at 187 yum which is beginning to show many of the 
features of the SCUBA 450 /tm map. Secondly, the dust opti¬ 
cal depth we derived from the Herschel data is quite similar to 
the 450 fim emission map. Finally, the luminosity image shows 
clearly that IRS 1 is the principal luminosity source in the re¬ 


Fig. 10 Temperatures of dust (red line) and gas along the HIFI cut. The 
gas temperature as derived fitting low-y lines (IRAM data, blue line) 
is systematically lower than the derived gas temperature fitting both 
low and high-/ lines (IRAM & HIFI data, black line). Tg^, (from both 
low and high-J lines) is systematically higher than T(i„j„ at least along 
this cut. The errors vary between ~5-15 % throughout the cut with the 
higher accuracy close to the central position. 


gion, followed in importance by IRS 2 and then IRS 3. There is 
no obvious luminosity peak within the area of strongest 450 fim 
emission, so this is probably a peak in the column densities of 
dust and gas. 

One uncertainty of the dust temperature may stem from the 
assumed spectral properties of the dust grains, in particular the 
assumed spectral index p. Recent observations of dust in the dif¬ 
fuse ISM with the Planck mission indicate a mean dust temper¬ 
ature of ~20 K and dust emissivity index of ~1.6 (Jones 2014). 
This is clearly higher than the value of P-\ assumed here, that 
is more appropriate for dense clouds (Ossenkopf & Henning 
1994). As our region is denser than the majority of those ob¬ 
served with Planck we can consider the value of P-l.b as an 
extreme upper limit for our case. When we fit a blackbody func¬ 
tion modified by a T ' ® to derive the dust temperatures we ob¬ 
tain temperatures that are lower by 5-20 K than the values from 
the adopted A '. We conclude that in our analysis we derive an 
upper limit of dust temperatures. 

4.3. Comparison of Dust and Gas temperatures. 

For a direct comparison with the gas kinetic temperature map, 
the dust temperature map was convolved to the same angular res¬ 
olution (~21"). The comparison between gas and dust tempera¬ 
tures from the low-T line analysis (Fig. 10 & Figure 13), shows 
a very similar spatial structure around the infrared sources. The 
higher-7 lines from HIFI though, reveal that the low-7 analysis 
underestimates the gas temperatures (see Sect. 4.1.3). Assuming 
that this trend applies to the entire region and not only along the 
cuts (Figure 10) we conclude that the gas temperature is higher 
than the dust in the entire region. 
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Fig. 11 Kinetic temperature (top) and column density (bottom) maps 
of S 140 (all lines > 3 RMS). The kinetic temperatures and column 
densities are higher towards the embedded infrared sources and lower 
throughout the rest of the cloud. The kinetic temperature rises again 
towards IF, but remains below the value towards the center. 
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The temperature difference of ~5-15 K revealed from the 
complete analysis lies outside the uncertainty range of the 2- 
7 K and indicates a more efficient gas heating even at densities 
> 10^'^ cm“^ where the two components are relatively well cou¬ 
pled. 

To illustrate the expected global temperature dependence, we 
show in Figure 14 the outcome of a simplified PDR model for 
S 140 computed with the KOSMA-t PDR code (Rollig et al. 
2006, 2013). It shows the temperature profile for a spherical 
PDR clump with a mass of 100 Mq, mimicking a plane-parallel 
PDR, and a surface density of 10^ cm“^ illuminated by an UV 
radiation field of 100 Draine fields (Draine 1978). For refer¬ 
ence, the green dash-dotted line shows the visual extinction Ay 
as a function of the depth into the cloud. The model shows that 
throughout almost all of the whole cloud, the gas temperature is 
higher than the dust temperature, by about 15 K for Ay = 0.1 and 
by only 2 K deep in the cloud. The small kink in the gas temper¬ 
ature around Ay = 1 stems from H 2 formation heating on PAHs 
and the temperature structure of the PAHs. Overall the model 
confirms that for a homogeneous medium, we expect the same 
temperatures for gas and dust within our observational measure¬ 
ment errors for visual extinctions into the cloud of more than 
about 0.5, i.e. for depths of one arcsecond and more at the dis¬ 
tance of S 140. The extended hot gas therefore cannot be ex¬ 
plained by a homogeneous extinction of the UV radiation but 
requires some clumpy structure allowing for a deeper UV pene¬ 
tration. 


Fig. 12 Left column: PACS/Spec continuum images (5x5 spatial array 
of 9" pixels regridded to a 1" grid) with contours of 37pm (SOFIA) 
emission overlaid showing the positions of IRS 1, 2, and 3. Right col¬ 
umn: top - total luminosity image with contours of fitted dust tempera¬ 
ture at 70, 65, 60, ...K; center - luminosity image with contours of fitted 
dust optical depth at levels of 0.1, 0.15, 0.2, and 0.25; bottom - 450 pm 
SCUBA image with same fitted dust optical depth contours overlaid. 

5. Effect of density gradients 

Assuming uniformly distributed densities and temperatures for 
each spatial position was a good approach for the extended gas 
and dust in S 140, but a detailed modeling is required for the 
more complicated dust and gas structures that are expected to be 
close to the infrared sources. The radiation we receive is a re¬ 
sult of various processes including the interaction with the sur¬ 
rounding dust. For a more precise physical approach, the radia¬ 
tive transfer calculations should include the original radiation, 
characteristics of the dust grains and the dust density distribu¬ 
tion. Density gradients have been revealed for S 140 in previous 
studies including Harvey et al. (1978), Giirtler et al. (1991), van 
der Tak et al. (2000), Maud & Hoare (2013). We first perform an 
advanced dust modeling taking into account dust density gradi¬ 
ents. Then we apply the best-fit derived dust model in order to 
predict the CO intensities using an advanced radiative transfer 
code (i.e. RATRAN; Hogerheijde & van der Tak 2000), as an at¬ 
tempt to test the accuracy of our dust model comparing with the 
gas observations and predictions. 
















Koumpia et al: Temperatures of Dust and Gas in S 140. 


63*2r00" 


63*20'00" 


63*19'00" 

o 

v 

o 


63*10'OO" 


63 * 17 ' 00 " 




z [pc] 


Fig. 14 Dust and gas temperatures in a spherical PDR clump with a 
surface density of 10^ cm“^ illuminated by an UV radiation field of 10^ 
Draine fields. Gas temperatures generally exceed the dust temperatures 
(see Rdllig et al. 2013, for a more general discussion). The green dash- 
dotted line shows the visual extinction Ay as a function of the depth into 
the cloud. 


Fig. 13 Gas kinetic temperature map towards IRAM field (low-T anal¬ 
ysis) with overplotted contours of dust temperatures (top) & zoomed in 
PACS field (bottom). The levels of contours are 55 K (inner) till 35 K 
(outer) using a step of 5 K. The colour bar indicates the values of ki¬ 
netic temperature. The 2 temperatures from this analysis seem to agree 
well but the gas kinetic temperatures are underestimated for ~7-12 K 
(Figure 10). 


5.1. Dust modeling - Approach 

It is clear from all the infrared and sub-mm images of the S 140 
cluster that there are a number of luminosity sources in the cen¬ 
tral arcminute and that the dust distribution is unlikely to be 
spherically symmetric. This implies that any simple model of the 
dust heating will be limited in its applicability. We have, how¬ 
ever, identified some goals to try to address in this study with 
simple radiative transfer models. First, since all the observational 
data suggest that IRS 1 is the dominant luminosity source, we 
believe that we can determine some rough properties of the dust 
distribution close to it by ignoring heating from IRS 2 and 3 and 
any of the other much lower luminosity objects that are part of 
the embedded stellar cluster around IRS 1. To this end, for the 
model comparison for IRS 1 we concentrate on two observables: 
(1) the SED of IRS 1, and (2) the spatial flux distribution from 
the center to the south of IRS 1 where there are no strong ob¬ 


vious nearby heating sources seen in the infrared that are likely 
to contribute significantly to the dust heating. A second goal is 
to try to understand whether the peaks in the 450 pm image can 
be due to some unusual dust distribution to the west and south¬ 
west of IRS 1 with no additional internal heat source or if some 
internal heating is required to explain them. For example, the re¬ 
cent observations and modeling by Maud & Hoare (2013) sug¬ 
gest an internal source within the strongest 450 pm peak to the 
southwest of IRS 1, i.e. SMM 1 (ARA: -14.22", ADec: -7.90", 
relative to IRS 1). 

The approach we took to addressing these goals involved 
several steps which we discuss in detail in the following sec¬ 
tions. First, we tried to find what properties of models were re¬ 
quired to provide a fit to the observed spectral energy distribu¬ 
tion (SED) from the central pixel centered on IRS 1 and to the 
source profile of IRS 1 to the south, assuming a single, spheri¬ 
cally symmetric dust distribution. This model was then used as a 
starting point for subsequent stages. The second step was to cre¬ 
ate a model made from hemispheres of two different spherically 
symmetric dust distributions. Though this model is simplified, it 
seems likely that far from the boundary region, such a combina¬ 
tion can give insight into the degree to which such a dust distri¬ 
bution might reproduce the observed source structure to the west 
and southwest of IRS 1 by assuming a higher column density in 
this region than to the south and east of IRS 1. For this compari¬ 
son we used the source profile extending southwest along a line 
from IRS 1 to and beyond the position of SMM 1 as given by 
Maud & Hoare (2013). As an alternative to the two-hemisphere 
model, we also attempted to fit the profiles along this line with 
two separate sources by using the superposition of two spheri¬ 
cally symmetric models with different central source luminosi¬ 
ties and very different dust distributions separated on the sky by 
~ 16" to attempt to reproduce the observed 1-D flux profile be¬ 
tween IRS 1 and the position of SMM 1. southwest of it. 

In testing these last two model variations, the two- 
hemisphere construction and the two-source construction, it be¬ 
came clear that the best results would likely result from a com¬ 
bination of both features. Such a model would also be most real¬ 
istic in light of the much higher dust column density to the west 
and southwest as indicated both by the SCUBA maps and our 
Herschel maps, together with the presence of the compact sub- 
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millimeter/millimeter source found by Maud & Hoare (2013). 
After constructing a few such models, we realized that with so 
many free parameters, it would be difficult to find the range 
of well-fitting parameters without testing the fits over a large, 
multi-dimensional grid. This process represented the final stage 
of modeling for IRS 1, and we now describe the process and 
results. Figure 15 shows a schematic diagram of the model for 
IRS 1 and the submm peak to the southwest. 


the range of values tested in a grid allowing all possible com¬ 
binations of these values (150,000 models). A preliminary run 
of ~ 800,000 models showed that the results were quite insen¬ 
sitive to the value of parameter 4, the radius where the slope of 
the density gradient changed in the high-column hemisphere, so 
it was fixed at a value of 300 x the inner radius. Similarly, the 
value of the density gradient in this inner region affected the final 
results only slightly, so this was fixed at a value 0.4 (p oc 
Therefore, for the final results there are seven free parameters. 


IRS 1 



Fig. 15 A schematic diagram in the place of the sky of the model for 
IRS 1 and the submm peak (ARA: -14.22", ADec: -7.90") to the south¬ 
west. The distance between the 2 positions is ~ 0.06 pc. 


For all the models we used the DUSTY code (Ivezic et al. 
1997) and converted the dimensionless output values to those 
appropriate for the assumed distance of S 140, 746 pc (Hirota 
et al. 2008) and the assumed luminosity (discussed below) for 
comparison with the observations. We tried models with two dif¬ 
ferent dust compositions; the first was a mixture of 90% DL sil¬ 
icates (Draine & Lee 1984) and 10% amorphous carbon, both of 
whose optical constants are distributed with DUSTY. The second 
was OH5 (Ossenkopf & Henning 1994) dust which a number 
of authors (e.g. van der Tak et al. 1999) have suggested may be 
more appropriate in dense star-forming regions where significant 
grain growth may have taken place. We found similar quality fits 
with comparable dust density distributions for both grain com¬ 
positions. The 1 yum optical depths were, however, of order twice 
as high with the DL/carbon dust than with the OH5 grain prop¬ 
erties due to the overall difference in the slope of the dust opti¬ 
cal depth between 1 and 100 fim. The output of the “observed” 
profile of each model at each wavelength was first interpolated 
onto a 3050x850 grid with spacing equivalent to 0.05". We then 
convolved this image with the appropriate observing beam to 
compare with the observed images (not those convolved to the 
187 yum resolution). Finally, we “observed” this image with a 
9" square beam that we moved across the image in 1" steps to 
compare with the observed images whose fluxes we also added 
into such a moving 9" square beam. For the shorter wavelength 
images we also applied the same 9" square pixel flux summation 
for comparison between the models and observations. 

5.2. Model Grid 

The model that we used is described by 9 parameters, four for 
the eastern half of IRS 1, four for the western half, and one, lu¬ 
minosity, for SMM 1. Table 2 lists these 9 parameters as well as 


Table 2 Model Parameters And Range for Grid Fitting With OH5 Dust. 


Parameter 

Value or Range 

Hot Source (IRSl) I.OxIO^Lq 

Density Gradient, r^'’ 

b^ 

d 

1 

o 

d 

Grain Size Distribution Slope 

3.5 

Grain Size Min/Max (ytim) 

O.l/l.O 

Dust Temp at Inner Radius 

1400K at 2.0x10*^ cm 


750 - 3000“ 

A., 

20 - 90 mag“ 

Cold Source (SMMl) 

Luminosity 

0 - 300Lo“ 

Photosphere Temperature 

2500 K 

Dust Temp at Inner Radius 

1200K 

Density Gradient, 

1.0 

^OKler/^inner 

3000 

Grain Size Distribution Slope 

3.5 

Grain Size Min/Max (ytim) 

O.l/l.O 

A., 

1000 mag 

Southwestern Profile 

Inner Density Gradient, r 

0.4 From R,„ to 300 Ry„ 

Outer Density Gradient, 

1.5 - -2.0 From 300R,„ to Ro„,“ 

^OKler/^inner 

750 - 3000“ 

Grain Size Distribution Slope 

3.5 

Grain Size Min/Max (ytim) 

O.l/l.O 

Dust Temp at Inner Radius 

1400K 

A., 

15-70 mag“ 


“Free parameter during modeling. 


In order to select the most likely models and to estimate the 
range of values that produce a reasonable fit, we characterized 
the model fits by the that we computed as follows. We pa¬ 
rameterized the observations and accompanying model fits into 
30 values of the SED and the source profile at various wave¬ 
lengths. Since the signal-to-noise ratio of almost all the observa¬ 
tions is quite high, the systematic uncertainties of the observa¬ 
tions (e.g. pointing) dominate the true uncertainties. We, there¬ 
fore, assigned somewhat arbitrary values to the assumed errors 
for these 30 values to drive the fitting to something “reasonably” 
close to the observations. Table 3 lists these properties of the 
observed data and uncertainties. 

As shown in the table, the parameters used to calculate the 
X^ of the fit included: the SED between 37 and 450 yum, the flux 
at the position of SMM 1, the relative flux at 37^50 jum at off¬ 
sets of 9" and 18"to either side of the position of IRS 1 (i.e. 
one pixel and two pixels), and finally, the observed luminosity at 
those four offset positions (not relative to the central pixel). We 
intentionally did not include the portion of the SED shortward 
of 37 fim, since the dust close to the central source of IRS 1 is 
likely to be distributed in a disk-like configuration which would 
lead to more flux escaping at shorter wavelengths than consis- 
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Table 3 Parameters and Assigned Relative Uncertainties Used to 
Compute of Fit To IRS 1. 


Parameter(s) 

Relative Uncertainty 

Peak Fy at 37, 450pm 

0.1 

PeakFy at 73, 125, 187pm 

0.05 

SMM 1 Fy at 450pm 

0.1 

Relative Fy at 37, 73, 125, 187, 450pm “ 

0.1 

Relative Luminosity “ 

0.1 


“The offsets were at -18", -9", 9"and 18"as appear 
along the dashed line in Fig 15. 


tent with the assumption of spherical symmetry. Likewise, we 
did not attempt to fit the total observed luminosity at the cen¬ 
tral pixel since much of the luminosity is emitted shortward of 
37 jjm. To estimate the range of model parameters that provided 
a “good” fit to the data, we assigned probabilities to each model 
equal to exp(-;(f^/2). 

The best-defined parameter values from the probability 
plots, are the mild density gradient away from IRS 1 of p oc r"® 
and the outer radius of 1500 times the inner radius. The optical 
depth is relatively well defined by ~ 40, though there is a 
better-defined joint probability that includes the optical depth 
and density gradient. The peak probability for the luminosity of 
SMM 1 is lOOLo, but there is a wide range of values with rea¬ 
sonable probability from ~ SOLq to BOOLq. The parameters of 
the dense material to the southwest of IRS 1 are not well defined 
for the most part. Interestingly the gradient of the material in 
the outer area (300 - 1500 x the inner radius) seems to have a 
most probable range similar to that found for the rest of the cloud 
around IRS 1, i.e. r~P with p ~Q- 0.5. There is a significant joint 
probability distribution between the outer gradient to the south¬ 
west and the luminosity of SMM 1. This is to be expected, since 
the effects of a steeper gradient can be largely compensated by 
increasing the luminosity of SMM 1. Although the most prob¬ 
able models have Ay = 40, the probability distribution suggests 
that the most probable range is 40 - 60, i.e. somewhat higher 
than the gradient to the south. 

Examples of the model fit profiles from one of the two low¬ 
est (2.7) fits as well as the SED are shown in Eigures 16-18. 
The model profiles shown in Eigure 16 display a less than per¬ 
fect fit at the longest wavelengths in our PACS observations, 140 
-190 pm. This effect has been driven by two parts of the fitting 
process, the fact that we did not include the relative flux at the 
central (IRS 1) position in the;^^^ and the fact that it proved quite 
difficult to find any models that dropped as slowly as the obser¬ 
vations at the ends of the profiles 18" away from the center. The 
properties of the models that do the best job fitting the profiles 
18" away from IRS 1 also drive the relative 187 pm profile to be 
too high at IRS 1 compared to its value at SMM 1. In some sim¬ 
ple tests we have found that if we added an artificial floor to the 
model fluxes beyond 125 pm, we could fit the longer wavelength 
profiles significantly better. 

5.3. IRS 2 and 3 

The most important factor affecting our attempts to model IRS 2 
and 3 is that the flux from our well-fitting models for IRS 1 is 
a significant fraction of the total flux at the positions of IRS 2 
and 3. This means that any models we make for IRS 2 and 3 will 
have an additional large uncertainty due to the uncertainty in the 
true contribution from IRS 1. Eor example, the model flux from 
IRS 1 at the position of IRS 3 is more than 50% of the total at 
wavelengths beyond 120 pm. We have therefore computed mod- 



-20 -10 0 10 -20 -10 0 10 
Offset from IRSI (Arcsec) Offset from IRSI (Arcsec) 

Fig. 16 Observed and model profiles at four PACS/Spec wavelengths 
normalized to the peak along a line from the south into the center of 
IRS 1 and then out along a line to the southwest extending to the SMM 1 
source. 



X (pm) 

Fig. 17 SED for a typical model fit with low;^^^. The solid line shows 
the total flux coming from the warm hemisphere, and the dotted line 
shows the SED for the fraction of the flux coming from the cold hemi¬ 
sphere. The plotted symbols show the observed and model fluxes for 
the fraction within the central 9" PACS pixel. 


els for roughly a dozen possible configurations for each source 
to develop some feeling for the range of likely parameters, but 
the uncertainties in our estimates are quite large. For IRS 3 we 
used a profile that begins on the east side of the source, and after 
passing through the center, goes straight to the south in order to 
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Fig. 18 Profiles along the same line used in Fig. 15 for the 450 /rm flux 
(top panel) and the total luminosity (bottom panel). Note the signiflcant 
mismatch in the central region which is due to not trying to fit the 5- 
30 yum luminosity since some of it likely arises from an axisymmetric 
disk rather than the spherical model used here. 



Fig. 19 Model and observed SED’s for IRS 2. The solid line shows the 
model fluxes for the hotter source to the west in the model, while the 
dashed line shows the fluxes for the colder source to the east. The sym¬ 
bols show the model and observed fluxes in the central 9" pixel as in 
Fig. 16. 


minimize the confusion from IRS 1. An additional complication 
is that IRS 2 is clearly elongated in the east-west direction in the 
mid-IR, and there is a clear position shift in the same direction 
between the mid-IR and our PACS data. This suggests that there 
are at least two relatively luminous objects heating the dust at 
the position of IRS 2 with different dust optical depths, higher 
to the east and lower to the west. The position shift between 12 
and 73 /rm is only a few arcsec, which is a fraction of our PACS 
native pixel size. Therefore, for IRS 2 we used a source profile 
from east to west and a model like that in the “two-hemisphere” 
modeling for IRS 1, with a hot, lower-optical-depth dust distri¬ 
bution on the western side and a cold, higher-optical-depth dis¬ 
tribution offset 4" to the east. Such a model will not accurately 
reproduce the source elongation, but should roughly fit the over¬ 
all SED and some part of the profiles. Figures 19-22 show the 
model SED and profile fits for some representative models that 
illustrate both the features that can be fit reasonably well as those 
that are difficult to fit. 

For both of these sources we did attempt to fit the mid- 
infrared fluxes in the SED, unlike for IRS 1. The very extended 
tail of emission on the west side of IRS 2 at 125 and 187 yum 
coincides well with the increase in optical depth seen in our data 
and the SCUBA maps in that area. We have not tried to fit this 
with these simple models. The bump seen on the eastern side of 
the long wavelength profiles for IRS 3 is coincident with the out¬ 
ermost pixel of our mapping, and we have likewise not attempted 
to fit it. 

The most well-defined results seem to be that: 1) the dust 
distribution around both sources is likely to be quite flat, and 2) 
the luminosities of the sources are likely in the range of 1000 - 
2 OOOL 0 , as shown in Table 4. The first could be consistent with 
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Relative Offset (Arcsec) Relative Offset (Arcsec) 

Fig. 20 Model and observed profiles for IRS 2 as in Fig. 15, but along a 
straight east-west line. 


very little dust around them that is associated with them, but 
rather mostly with the extended distribution around IRS 1. 

5.4. Gas modeling 

In order to test our best fit DUSTY model we ran the Monte 
Carlo radiative transfer code RATRAN (Hogerheijde & van der 
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Fig. 21 Model and observed SEDs for IRS 3. 
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Fig. 22 Model and observed profiles for IRS 3 along a line extending 
into the source center from the east, and then outward to the south. 


Tak 2000) which treats the physical structure of the sources in¬ 
cluding temperature and density gradients. RATRAN estimates 
the local radiation field at all line frequencies taking into ac¬ 
count the radiation field from every other position in the cloud. 
We ran RATRAN applying the best fit DUSTY model and as¬ 
suming the same temperature for dust and gas towards IRS 1 in 
order to model the CO 1-0, CO 2-1, CO 9-8 and isotopologue 
lines and compare them to the outflow-subtracted component 
from the observations (Figure 23). The outflow emission was 
subtracted since the DUSTY models focus on the bulk of the 
quiescent dense gas where the outflow cavities are negligible. 


Table 4 Model Parameters For Symmetric Fit With OHS Dust to IRS 3 
and “Two-Hemisphere” Fit To IRS 2. 


Parameter 

IRS 2 

IRS 3 

Single Gradient 

Eastern Side 

Symmetric 

Source Luminosity 

2000 Lo 

1300 Lo 

Dust Temp at Inner Radius 

1400K 

1400K 

Density Gradient, 

0.0 

0.0 

Grain Size Distribution Slope 

3.5 

3.5 

Grain Size Min/Max (pm) 

O.l/l.O 

O.l/l.O 

^ouierf^inner 

3000 

1500 

A. 

30 mag 

20 mag 

Double Gradient 

Western Side 


Inner Density Gradient, r^’’ 

0.3 From R,>, to 300 R,„ 


Outer Density Gradient, r^’’ 

0.3 From 300R,„ to Ro„, 


^ouier/^inner 

1500 


A„ 

110 mag 



For the outflow subtraction we applied a 2 component Gaussian 
fit to the observed lines towards IRS 1 followed by the subtrac¬ 
tion of the broad component. 

For our calculations we defined a grid of 18 spherical shells 
for the east (lower density) hemisphere (constant p oc "^) and 
22 for the west (high density) hemisphere applying the density 
gradients as obtained from the DUSTY model. The inner and 
outer radius were taken from the DUSTY model and were set to 
2.08x10'“* cm and 3.12x10'^ cm respectively, while the inner H 2 
density was set to 9x10^ cm“^. This value was derived from the 
best fit dust model assuming that the gas is entirely molecular 
and using a mean gas mass per hydrogen of 1.4 amu and a gas- 
to-dust ratio of 100. 

Dust continuum radiation is taken into account using the 
same OHS opacities as in the DUSTY model. We assume a static 
envelope without infall or expansion, a fixed CO abundance of 
10 * and a turbulent line width of 3.5 km s *. 

Figure 23 shows the lines of CO 1-0, CO 2-1 and CO 
9-8 and isotopologues as observed towards IRS 1 overplotted 
with the convolved (~21") synthetic emission as calculated with 
RATRAN. The red color represents the resulting line profiles 
when using the lower density east hemisphere from the DUSTY 
models, while the blue color represents the high density west 
hemisphere. With the exception of the low-7 transitions of the 
isotopologues, we observe no significant differences between the 
two cases. The fit between observed and modeled lines shows 
that the low-7 lines are reproduced but the higher-7 lines are 
significantly underestimated, especially for the isotopologues. 
The model cannot reproduce the highly excited lines tracing high 
column densities of warm and very dense gas. An adapted kine¬ 
matic structure and/or dumpiness can potentially remove the 
modeled self-absorption dip of the optically thick *^CO lines 
but treating these effects are beyond the scope of this work. 
Models with Tg^j > T^ust, as inferred from the simple analysis in 
Sect. 4.3, would probably be able to reproduce the high-7 lines 
without significantly affecting the low-7 lines. Local density en¬ 
hancements can be an alternative solution, if the total mass of 
the cloud is conserved, such as in a clumpy or disk/outflow ge¬ 
ometry. 
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Fig. 23 Outflow subtracted spectra of CO 1-0, CO 2-1, CO 9-8, '^CO 
1-0, ‘^CO 10-9, C‘*0 1-0 and C‘'*0 9-8 (black) as observed to¬ 
wards IRS 1, overplotted with the synthetic emission as calculated 
with RATRAN having applied the lower density east hemisphere from 
DUSTY models (red) and the high density west hemisphere (blue). We 
observe no significant differences between the two cases with the excep¬ 
tion of the low-J transitions of the isotopologues. While the model fits 
well the observations for the lower transition lines, it fails to reproduce 
the observed higher transition lines. 


6. Discussion and conciusions 

Based on PACS, HIFI and IRAM data, we find a large range 
of dust and gas temperatures in the S 140 region. The warmest 
gas (~ 55 K) is around the three infrared sources, while the sur¬ 
rounding environment is characterized by colder temperatures 
(~25-35 K). The gas temperature increases towards the south¬ 
west edge of the cloud (IF) reaching values of ~40^5 K. This 
rise of temperature is the result of the UV heating from the ex¬ 
ternal BOV star HD211880 that lies ~T southwest of the edge 
of the cloud. Unfortunately, we cannot compare the gas and the 
dust temperatures towards IF, because of the smaller area that 
our continuum dust observations cover. 

We find that the dust density gradient around IRS 1 is likely 
to be shallower than the best fits found in some earlier models in 
order to explain the spatial structure we find in the far-infrared. 
The strong emission to the west and southwest of IRS 1 at wave¬ 
lengths longward of 100 jim can possibly be explained solely by 
heating from IRS 1 together with a strongly increasing density to 
the west and southwest. We do, however, find a significantly bet¬ 
ter fit by interpreting the emission as arising from two or more 
separate sources with very different amounts of dust around the 
sub-mm source compared to IRS 1. The luminosity of the in¬ 
ternal source at the sub-mm peak is likely to be few x 10 ^ Lq 
in this model. This model with an internal luminosity source at 
the sub-mm peak is consistent with the study of Maud & Hoare 
(2013) who find that the two strongest 1.3 and 2.7 mm emission 


peaks are the position of IRS 1 and SMM 1. Since our models 
did not include outflow cavities or dumpiness, it is possible that 
this conclusion should be revisited after more extensive model¬ 
ing. 

The gas temperature analysis including both high and low-7 
lines was possible only along the HIFI cut and revealed a system¬ 
atic excess of gas temperatures against dust temperatures. This 
result indicates a more efficient gas heating even at densities > 
10 "^ ^ cm“^ where the gas and dust are usually expected to be 
in thermal equilibrium. New SOFIA/GREAT mapping observa¬ 
tions of CO 13-12 and 16-15 confirm this excess increasing the 
gap between the two even more (Ossenkopf et al. submitted). 
The high-7 lines show that a limitation to the low-7 analysis 
(IRAM) causes an underestimate of the gas temperatures by ~1- 
12 K and thus provides only a lower limit of the temperatures in 
the cloud. Thus it is most likely that the gas temperature is higher 
than the dust also in the whole field. Unfortunately, the higher-7 
lines do not cover the entire region so that we cannot prove this. 
The detailed gas modeling also indicates that DUSTY results are 
good estimates of the temperature and density structure around 
IRS 1 for the colder gas but there is an excess of warm dense gas 
that cannot be reproduced by these models (i.e. CO 9-8). The at¬ 
tempt to fit those higher transition lines by modeling hotter gas 
than dust and/or applying different geometries (i.e. clumps, disk 
geometry) that would increase the density locally but not the to¬ 
tal amount of mass is in our future plans. 

The observed gas temperature excess cannot be explained 
by a higher cosmic ray ionization. Van der Tak & van Dishoeck 
( 2000 ) reported a cosmic ray ionization rate of < 10 “*® s“*, 
which is too low to provide a major gas heating in S 140. Some 
gas heating may be due to outflow activities. In the RADEX and 
RATRAN analysis, we focused on quiescent gas and limited the 
contribution from protostellar outflows. However, it is possible 
to have high velocity motions perpendicular to the line of sight. 
In addition, the larger difference between gas and dust appears 
in points close to the observed outflow by Preibisch & Smith 
(2002). A collection of shocks, longitudinal and transversal to 
the walls of the possible inner cavity around the IR sources will 
contribute to the emission at the ambient velocity. Oblique shock 
components are thus a plausible scenario for a somewhat more 
efficient heating of the gas compared to the dust. 

The most probable scenario appears to be the deep UV radi¬ 
ation in a clumpy medium. Dust is heated by UV and IR while 
gas is heated by processes driven only by UV with the main one 
being the photo-electric heating (Hollenbach & Tielens 1997; 
Rollig et al. 2013). In the vicinity of any stellar or protostellar 
source the gas is much hotter than the dust when dust and gas 
are not coupled. Dust-gas collisions try to equilibrate them if 
the density is high enough. The dust quickly attenuates the UV 
so that the gas becomes colder when going away from IRS 1. 
The penetration of the IR is deeper (lower extinction in IR com¬ 
pared to UV), so that the decay of the temperature of the dust is 
shallower than that of the gas. However, as we observe that the 
gas remains hotter than the dust over a relatively large distance 
from IRS 1, the UV seems to penetrate much deeper into the 
cloud than expected from a homogeneous medium. The natural 
explanation is a clumpy medium where there are always enough 
rays between the clumps that still allow UV propagation to keep 
the gas warm. 

Previous studies already report S 140 to be clumpy (e.g. 
Kramer et al. 1998). Battersby et al. (2014) performed a similar 
study in order to determine the relationship between gas and dust 
in a massive star-forming region (G 32.02-1-0.05) by compar¬ 
ing the physical properties derived from each, using NH 3 tran- 
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sitions. In that study similar temperature differences have been 
reported with the dust temperatures being lower than gas tem¬ 
peratures (by a few K) in the quiescent region, indicating that 
gas and dust might be not well-coupled in such environments. 
Differences between gas and dust temperatures in several envi¬ 
ronments with high densities have been reported in more studies 
(i.e. Papadopoulos et al. 2011; Hollenbach & Tielens 1999). 

Our observations trace warm and cold dust but also warm 
and less warm gas since we have multiple transitions of CO. The 
low-J lines tend to show the surface of clumps since they are op¬ 
tically thick but the optically thin isotopologues '^CO 1-0 and 
C'^O 1-0 carry information from deeper positions. Furthermore 
the higher transitions CO 9-8 and '^CO 10-9 trace warmer gas 
but they also arise from deeper positions. The fact that we have 
optically thin lines and the dust emission is optically thin con¬ 
firms that we trace the whole line of sight both in gas and dust 
and thus the observational limitations can probably be excluded. 
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